Instructions to access sf object mapping from Chris’ GitHub: https://github.com/cfree14/wcfish
# Run if you don't already have devtools installed
# install.packages("devtools")
# Run once devtools is successfully installed
# devtools::install_github("cfree14/wcfish", force=T)
# library(wcfish)
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ tibble 3.1.8 ✔ dplyr 1.1.0
## ✔ tidyr 1.2.1 ✔ stringr 1.4.1
## ✔ readr 2.1.3 ✔ forcats 0.5.2
## ✔ purrr 0.3.5
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks plotly::filter(), stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(here)
## here() starts at /Users/clairegonzales/Documents/R projects/chap_2_sitingstudy/chap2
library(sf)
## Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE
#library(tmap)
# library(paletter)
Getting the blocks shapefile from Chris’ Github:
# See dataset metadata by typing a question mark and the dataset name
# ?wcfish::blocks
# Store the dataset in an object of your own choosing as follows:
blocks_sf <- wcfish::blocks %>%
drop_na(block_lat_dd, block_long_dd)
# %>%
# st_as_sf(coords = c("block_long_dd", "block_lat_dd")) #take lat and long and convert to spatial points
################################################################################
# FOR ACTUAL MAPPING
# Use the ESRI .shp file from: https://chrismfree.com/california-fisheries-data/ca-fishing-blocks/
blocks_shp <- read_sf(here("data_fromChris", "CA_OR_WA_comm_fishing_blocks", "CA_OR_WA_comm_fishing_blocks.shp"))
Read in catch data from Chris
catchdata <- read_csv(here("data_fromChris", "gonzalez_data.csv"))
## New names:
## Rows: 640 Columns: 8
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (1): confidential dbl (7): ...1, block_id, nvessels, nfishers, nbusinesses,
## value_usd, landing...
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...1`
Now that we have both the catch data from Chris and the blocks shapefile, I am going to run these through a spatial function
# looking at blocks spatially
ggplot(data = blocks_sf, aes(x = block_long_dd, y = block_lat_dd)) +
geom_point(shape = 15, fill = NA, aes(size = block_sqkm)) +
theme_minimal()
Making the catch data ready for spatial mapping
catchdata_sf_merge <- merge(blocks_sf, catchdata) %>%
filter(block_state == "California",
block_type %in% c("Inshore", "Midshore")) %>%
mutate() %>%
select(block_id, block_type, block_sqkm, block_long_dd, block_lat_dd, nvessels:geometry) %>%
drop_na(block_long_dd, block_lat_dd) %>%
st_as_sf(coords = c("block_long_dd", "block_lat_dd")) #take lat and long and convert to spatial points
Coming back later to add a column of the z transformed data in R. This seems easier to do here than in ArcPro.
a <- min(catchdata_sf_merge$landings_lb, na.rm = TRUE)
b <- max(catchdata_sf_merge$landings_lb, na.rm = TRUE) + 1
m <- (a + b)/2
# median(catchdata_sf_merge$landings_lb, na.rm = TRUE)
catchdata_sf_merge_z <- catchdata_sf_merge %>%
mutate(z_lndngs_ = case_when(
landings_lb <= a ~ 1,
(landings_lb > a & landings_lb <= m) ~ (1-(2*(((landings_lb-a)/(b-a))**2))),
(landings_lb > m & landings_lb < b) ~ (2*(((landings_lb - b)/(b-a))**2)),
landings_lb >= b ~ 0,
TRUE ~ as.numeric(NA)
))
# trying to normalize the data by taking the natural log
catchdata_sf_log <- catchdata_sf_merge %>%
mutate(ln_landngs = log(landings_lb))
# WORKED. But also now everything is negative soooo let's see how that looks
Trying a linear transformation of the data
slope <- (1/-(b-a))
slope2 <- (1/-b)
y_int <- 1
catchdata_sf_merge_linear <- catchdata_sf_merge %>%
mutate(linear_landings = (landings_lb*slope2)+1)
## See bottom chunk for plots and plotly of this transformed data
Option 4 - log transformation of data and then Z transformation
a_ln <- min(catchdata_sf_log$ln_landngs, na.rm = TRUE)
b_ln <- max(catchdata_sf_log$ln_landngs, na.rm = TRUE) + 1
m_ln <- (a_ln + b_ln)/2
catchdata_sf_merge_ln_z <- catchdata_sf_log %>%
mutate(z_lndngs_ln = case_when(
ln_landngs <= a_ln ~ 1,
(ln_landngs > a_ln & ln_landngs <= m_ln) ~ (1-(2*(((ln_landngs-a_ln)/(b_ln-a_ln))**2))),
(ln_landngs > m_ln & ln_landngs < b_ln) ~ (2*(((ln_landngs - b_ln)/(b_ln-a_ln))**2)),
ln_landngs >= b_ln ~ 0,
TRUE ~ as.numeric(NA)
))
Let’s take a look at the transformed data
#scatter plot of ln and z transformed data, against original catch data
ggplot(data = catchdata_sf_log, aes(x = landings_lb, y = ln_landngs)) +
geom_point() +
theme_minimal()
## Warning: Removed 1 rows containing missing values (`geom_point()`).
# histogram of original catch data
ggplot(data = catchdata_sf_log, aes(x = landings_lb)) +
geom_histogram() +
theme_minimal() +
xlim(0, 2500000)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 27 rows containing non-finite values (`stat_bin()`).
## Warning: Removed 2 rows containing missing values (`geom_bar()`).
# histogram of ln and z transformed data
ggplot(data = catchdata_sf_log, aes(x = ln_landngs)) +
geom_histogram() +
theme_minimal()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 1 rows containing non-finite values (`stat_bin()`).
# we need to tell r what the coord reference system is using `st_crs`
#st_crs(catchdata_sf_merge) <- 4326
#UPDATE
# Doing this with the transformed data instead
st_crs(catchdata_sf_merge_z) <- 4326
Add blocks shapefile to map
#
# blocks_shp_transform <- st_transform(blocks_shp, 4326)
#
# ggplot(data = blocks_shp_transform) +
# geom_sf()+
# theme_minimal()
Combine catchdata and block shapfile map
# map <- ggplot() +
# geom_sf(data = blocks_shp_transform) +
# geom_sf(data = catchdata_sf_merge, na.rm = TRUE,
# aes(color = landings_lb, size = block_sqkm)) +
# theme_void()
#trying this a different way
catchdata_sf_merge_transform <- st_transform(catchdata_sf_merge_z, 4326)
#map of landings
ggplot(data = catchdata_sf_merge_transform) +
geom_sf(aes(fill = landings_lb)) +
guides(color = guide_legend(reverse = TRUE)) +
theme_minimal()
#map of value
ggplot(data = catchdata_sf_merge_transform) +
geom_sf(aes(fill = value_usd))+
theme_minimal()
# looking at data
ggplot(data = catchdata_sf_merge_transform, aes(y = landings_lb)) +
geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 1 rows containing non-finite values (`stat_bin()`).
# Looks like the data is strongly skewed left, which is why the legend is so not useful
Write the catch data as a shapefile so that we can upload to ArcPro
st_write(catchdata_sf_merge_transform, here("data", "catchdata_output", "catchdata_output.shp"), append = FALSE)
## Warning in abbreviate_shapefile_names(obj): Field names abbreviated for ESRI
## Shapefile driver
## Deleting layer `catchdata_output' using driver `ESRI Shapefile'
## Writing layer `catchdata_output' to data source
## `/Users/clairegonzales/Documents/R projects/chap_2_sitingstudy/chap2/data/catchdata_output/catchdata_output.shp' using driver `ESRI Shapefile'
## Writing 547 features with 12 fields and geometry type Multi Polygon.
OPTIONS: Data that is transformed to Z-shaped membership function OR Linear transformation
plotz <- ggplot(data = catchdata_sf_merge_z, aes(x = landings_lb, y = z_lndngs_)) +
geom_point() +
theme_minimal()
plotz_plotly <- ggplotly(plotz)
plotz_plotly
#scatter plot of ln and z transformed data, against original catch data
plot1 <- ggplot(data = catchdata_sf_log, aes(x = landings_lb, y = ln_landngs)) +
geom_point() +
theme_minimal()
plot1_plotly <- ggplotly(plot1)
plot1_plotly
plot2 <- ggplot(data = catchdata_sf_merge_linear, aes(x = landings_lb, y = linear_landings)) +
geom_point() +
theme_minimal()
plot2_plotly <- ggplotly(plot2)
plot2_plotly
plot4 <- ggplot(data = catchdata_sf_merge_ln_z, aes(x = ln_landngs, y = z_lndngs_ln)) +
geom_point() +
theme_minimal()
plot4_plotly <- ggplotly(plot4)
plot4_plotly